#include <stdio.h>
#include <stdlib.h>

const int threadsPerBlock = 32;

typedef struct
{
    int J_ORDER;            // Ordem da matriz
    int J_ROW_TEST;         // Linha para ser avaliada
    float J_ERROR;          // Erro permitido
    int J_ITE_MAX;          // Número máximo de iterações
    float *MA;              // Matriz A do sistema         
    float *MB;              // Matriz B do sistema

    float *MA_J_ROW_TEST;   // Armazena a linha a ser avaliada
    float MB_J_ROW_TEST;    // Armazena o resultado da linha a ser avaliada
}Jacobi;


/*
*   Aloca as matrizes da struct 'Jacobi' e inicializa seus valores pela leitura do arquivo.
*/
void inicia_Jacobi( FILE *pFile, Jacobi *J )
{
    int i, j;
    int ORDER;
    int LINHA;

    fscanf( pFile, "%d", &(J->J_ORDER) );
    fscanf( pFile, "%d", &(J->J_ROW_TEST) );
    fscanf( pFile, "%f", &(J->J_ERROR) );
    fscanf( pFile, "%d", &(J->J_ITE_MAX) );

    ORDER = J->J_ORDER;

    J->MA = (float *) malloc( (ORDER * ORDER) * sizeof(float) );
    J->MB = (float *) malloc( ORDER * sizeof(float) );

    for(i=0; i<ORDER; i++)
    {
        LINHA = i*ORDER;

        for(j=0; j<ORDER; j++)
            fscanf( pFile, "%f", &(J->MA[LINHA+j]) );
    }

    for(i=0;i<ORDER;i++)
        fscanf( pFile, "%f", &(J->MB[i]) );

    J->MA_J_ROW_TEST = (float *) malloc( ORDER * sizeof(float) );
}

/*
*   Prepara as matrizes MA e MB para a aplicação do método de Jacobi-Richardson.
*
*   Antes de manipular a linha das matrizes, verificamos se ela corresponde a linha que será feita o teste no final.
*   Caso seja, armazenamos seus valores.
*
*   Cada linha das matrizes, então, é dividida pelos respectivos valores da diagonal da mesma linha.  
*
*   No vetor 'x' é armazenado os valores iniciais do método Jacobi-Richardson. 
*/
void prepara_Jacobi(Jacobi *J, float *x, int id)
{
    int i, j;
    int LINHA;
    int ORDER;
    float div;

    ORDER = J->J_ORDER;

    for(i=0; i<ORDER; i++)
    {
        LINHA = i*ORDER;

        /* Armazena a linha caso seja a linha do teste final (J_ROW_TEST) */
        if(i == J->J_ROW_TEST)
        {
            for(j=0; j<ORDER; j++)
            {
                J->MA_J_ROW_TEST[j] = J->MA[LINHA+j];
            }

            J->MB_J_ROW_TEST = J->MB[i];
        }

        /* Armazena a diagonal da linha */
        div = J->MA[LINHA+i];

        /* Divisão da linha pela diagonal */
        for(j=0; j<ORDER; j++)
        {
            J->MA[LINHA+j] /= div;
        }

        J->MB[i] /= div;

        /* Zera a diagonal da linha pois é ignorada no cálculo */
        J->MA[LINHA+i] = 0;

        /* Armazena no vetor 'x' os valores iniciais do método Jacobi-Richardson */
        x[(id*ORDER)+i] = J->MB[i];
    }
}

/*
*   @Método Jacobi-Richardson
*
*   A partir dos valores anteriores armazenados em x[id][...], realizamos uma nova iteração
*   e armazenamos os novos valores novamente em x[1][...] ou x[0][...], dependendo do valor de 'id'.
*
*   O erro a partir dos valores de 'x' é então retornado.  
*/

__global__ void cuda_jacobiRichardson(float *MA, float *MB, float *x, int *id)
{
    int i;

    __shared__ float cache[threadsPerBlock];
    int id_MA = threadIdx.x + blockIdx.x * gridDim.x;
    int id_MB = blockIdx.x;
    int id_x = threadIdx.x;
    int cacheIndex = threadIdx.x;
    int index = *id;
    float temp = 0;

    int N = blockDim.x;

    while(id_x < gridDim.x)
    {
        temp -= MA[id_MA] * x[id_x + gridDim.x * index]; //improve

        id_MA += N;
        id_x += N;
    }

    cache[cacheIndex] = temp;

    __syncthreads();

    if(index == 0) 
        index = 1;
    else 
        index = 0;


    i = N/2;
    while(i != 0)
    {
        if(cacheIndex < i)
        {
            cache[cacheIndex] += cache[cacheIndex + i];
            __syncthreads();
        }

        i /= 2;
    }

    if(cacheIndex == 0)
    {
        x[blockIdx.x + gridDim.x * (index)] = cache[0] + MB[id_MB];
    }

    if(blockIdx.x == 0 && threadIdx.x == 0)
    {
        if(*id == 0)
            *id = 1;
        else
            *id = 0;
    }
}

float error(float *x, int id, int ORDER)
{
    int i;
    double aux;
    double erro1 = 0;
    double erro2 = 0;

    for(i=0; i<ORDER; i++)
    {
        if(id == 1)
        {
            /* Verifica o erro */
            aux = x[ORDER+i] - x[i];

            if(aux < 0) aux *= -1;

            if(aux > erro1 ) erro1 = aux;

            aux = x[ORDER+i];

            if(aux < 0) aux *= -1;

            if(aux > erro2 ) erro2 = aux;
        }
        else
        {
            /* Verifica o erro */
            aux = x[i] - x[ORDER+i];

            if(aux < 0) aux *= -1;

            if(aux > erro1 ) erro1 = aux;

            aux = x[i];

            if(aux < 0) aux *= -1;

            if(aux > erro2 ) erro2 = aux;
        }
    }

    /* Retorna o erro */
    return erro1/erro2;
}

/*
*   Imprime a saída pelo padrão do projeto.
*/
void saida(Jacobi J, float *x, int id, int iter)
{
    int i;
    int LINHA;
    float aux=0;

    LINHA = id * J.J_ORDER;
    
    for(i=0;i<J.J_ORDER;i++)
    {
        aux += J.MA_J_ROW_TEST[i] * x[LINHA+i];
    }

    printf("Iterations = %d\n", iter);
    printf("RowTest: %d => [%f] =? %f\n", J.J_ROW_TEST, aux, J.MB_J_ROW_TEST);  
}

int main()
{
    int iter = 0;   // Número de iterações
    int *id;        // Índice para controlar os resultados anteriores e novos da matriz 'x' (Será sempre 0 ou 1) ´
    float *x;       // 'x' irá guardar os novos valores e os anteriores de cada iteração do método numérico 
    float erro;     // Armazena o erro do método
    FILE *pFile;
    Jacobi J;

    /* Arrays que serão armazenadas no device */
    int *dev_id;
    float *dev_MA, *dev_MB;
    float *dev_x;

    pFile = fopen ("matriz1000.txt","r");

    /* Inicializa a struct */
    inicia_Jacobi(pFile, &J);

    /* Aloca 'x' onde será armazenado os resultados das iterações */
    x = (float *) malloc( (J.J_ORDER * 2) * sizeof(float) );

    /* Aloca 'id' que controla o vetor 'x' ao armazenar os novos valores */
    id = (int *) malloc( sizeof(int) );

    *id = 0;

    /* Prepara a struct */
    prepara_Jacobi(&J, x, *id);

    /* Aloca memória na GPU */
    cudaMalloc( (void**)&dev_MA, (J.J_ORDER * J.J_ORDER) * sizeof(float) );
    cudaMalloc( (void**)&dev_MB, (J.J_ORDER) * sizeof(float) );
    cudaMalloc( (void**)&dev_x, (2 * J.J_ORDER) * sizeof(float) );
    cudaMalloc( (void**)&dev_id, sizeof(int) );

    /* Copia as arrays MA, MB, x e id para a GPU */
    cudaMemcpy( dev_MA, J.MA, (J.J_ORDER * J.J_ORDER) * sizeof(float), cudaMemcpyHostToDevice );
    cudaMemcpy( dev_MB, J.MB, (J.J_ORDER) * sizeof(float), cudaMemcpyHostToDevice );
    cudaMemcpy( dev_x, x, (2 * J.J_ORDER) * sizeof(float), cudaMemcpyHostToDevice );
    cudaMemcpy( dev_id, id, sizeof(int), cudaMemcpyHostToDevice );

    /* @Jacobi-Richardson */
    do
    {
        /* Cada bloco irá tratar uma linha do sistema */
        cuda_jacobiRichardson<<<J.J_ORDER,threadsPerBlock>>>( dev_MA, dev_MB, dev_x, dev_id );

        /* Copia os novos valores de x para a CPU e verifica o erro */
        cudaMemcpy( x, dev_x, (2 * J.J_ORDER) * sizeof(float), cudaMemcpyDeviceToHost );
        cudaMemcpy( id, dev_id, sizeof(int), cudaMemcpyDeviceToHost );

        iter++;

        erro =  error( x, *id , J.J_ORDER );

    /* Sai do método caso o erro seja inferior a J.J_ERROR ou ultrapasse o número máximo de iterações permitidas*/
    }while(erro > J.J_ERROR && iter < J.J_ITE_MAX);

    cudaMemcpy( x, dev_x, (2 * J.J_ORDER) * sizeof(float), cudaMemcpyDeviceToHost );
    cudaMemcpy( id, dev_id, sizeof(int), cudaMemcpyDeviceToHost );

    saida(J, x, *id, iter);

    printf("\n\nFIM");

    cudaFree( dev_MA );
    cudaFree( dev_MB );
    cudaFree( dev_x );
    cudaFree( dev_id );

    return 0;
}
